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Simulations of the early evolution of induced unstable disturbances in a Mach 6 hyper- 
sonic boundary layer are presented for the purposes of validating a high-order discontinuous 
Galerkin Navier-Stokes simulation code. The simulations are modeled after experiments 
performed in the Boeing/AFOSR Mach-6 Quiet Tunnel at Purdue University. In these 
experiments, a forcing mechanism was employed to induce reproducible disturbances in a 
hypersonic boundary layer providing for the controlled study of the growth and breakdown 
of these disturbances into turbulent spots. Simulations revealed that the form and strength 
of the excitation can greatly influence the growth of the disturbance. In particular, at large 
forcing amplitude, the simulated forcing produces large advecting transients that appear to 
enhance the growth of the wave packet, relative to that of low amplitude forcing. Simula- 
tion results with large amplitude forcing agree well with experimental results while results 
from low amplitude forcing agree with linear stability theory. 


I. Introduction 

Experiments recently performed in the Boeing/AFOSR Mach-6 Quiet Tunnel (BAM6QT) at Purdue 
University 1-3 studied the growth and breakdown of disturbances in a well controlled environment. This 
experiment was chosen as a validation case for a high-order three-dimensional Navier-Stokes simulation code 
based on the discontinuous Galerkin (DG) methodology. Prior work validating a two-dimensional variant of 
this code demonstrated super-convergence of unsteady viscous flow in the wake of a cylinder. 4 Some recent 
work to improve robustness of the DG method in the presence of shocks is described in Ref. 5. That work 
also describes some preliminary simulations that employed an axisymmetric variant of the current code. 

In the experiment, illustrated in Fig. 1 taken from Ref. 1, an electrode induces a disturbance that 
grows quickly into a nonlinear turbulent spot. Pressure is measured on the tunnel wall at several points 
downstream of the electrode. Typical ensemble-averaged results from the Purdue experiment are shown in 
Fig. 1(c). When operated in quiet flow mode, BAM6QT produces a relatively thick laminar boundary layer 
with a long transition region that enables accurate and detailed measurements with the available pressure 
sensors. Ideally, a study of naturally occurring disturbances would be preferred. However, the naturally 
occurring disturbances form in random locations relative to the pressure sensors and vary in strength making 
meaningful measurement difficult. The experiment employs an electrode forcing mechanism resulting in the 
production of consistent and repeatable wave packets. 

In early experimental work 1,2 the electrode was thought to produce a spark, presumably of very short 
duration; and the electrode was referred to as a “spark perturber.” Guided solely by pressure measurements 
just downstream of the electrode and the assumption of a forcing of short duration, early axis-symmetric 
simulations 5 attempted to reverse engineer a numerical forcing that would produce a similar downstream 
disturbance. However, the short duration of the forcing required an extremely large amplitude to produce any 
measurable disturbance. The large amplitude forcing posed a challenge to the robustness of the numerical 
simulation method, which would fail during the initial transient phase if the amplitude was too large. Efforts 
to improve the robustness of the method in the presence of shock waves 5 were moderately successful, allowing 
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(b) Schematic of experimental setup with perturber and sensor locations denoted on the x-axis. 



Figure 1. Experiment configuration and typical results from Ref. 1 
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the use of larger amplitudes than otherwise possible. Recent experimental work 3 measured characteristics of 
the forcing and determined that it was of much longer duration and did not normally involve a spark at all. 
Instead, the electrode is described as a “pulsed glow perturber.” The longer duration reduces the required 
amplitude of the forcing and results in a much more benign initial transient. 

The next section of this paper describes the governing equations and the DG discretization method used 
to solve them. This is followed by a description of the computational domain and the forcing function used 
to model the pulsed glow perturber. Finally, results of three-dimensional simulations are presented and 
compared to both experimental results and to linear stability analysis. 

II. Governing Equations and DG Formulation 

The nondimensional compressible Navier-Stokes equations for a perfect gas in conservation form, using 
tensor index notation are: 


dp djpuj) = 

dt dxj 

dpe d ( huj - UiT Zt j + qj) _ 

dt dxj 

dpui d(pUiUj + 5ijP - Ti tj ) 

dt dxj 


( 1 ) 

( 2 ) 

i = 1,2,3 (3) 


where p is the density, e is the internal energy per unit mass, P is the pressure, iq is the component of 
velocity in the Cartesian coordinate direction Xi, and Sij is the Kronecker delta. The quantities h, qi , and 
Tij are the enthalpy, heat flux, and shear stress terms, respectively; and are given by: 


h = pe + P, 


7 1 p dT 

^ 7 — 1 Pr Re r dxi ’ 

p ( dui duj 2 du n \ 

Re r V dxj dxj l ’ j 3 dx n ) ’ 

where Pr is the Prandtl number, and Re r is the Reynolds number based on the reference state of the 
nondimensionalization, and T is the temperature given by T = P/ p = (7 — l)(e — upiij 2). The length and 
the thermodynamic variables have been nondimensionalized with respect to a prescribed reference state: 
Xi = Xi/L r , p = p/p r , P = P/Pr 1 T = T/T r , and p = p/p r where ( A ) denotes dimensional quantities, 
and the subscript r denotes the reference state. Other variables are normalized with respect to derived 

reference states as follows: u = u/u r , t = t/(L r /u r ), e = e r /u and Re r = p r u r L r / p r where u r = \j Pr/ Pr- 

Sutherland’s formula is used to evaluate both p r as a function of T r and p as a function of T. 

DG is applied to each of Eqs. 1 3 in essentially the same manner, however, the inviscid and viscous terms 
are treated differently. To facilitate the following discussion, each equation is cast in the general form: 

^+V-(F < -F„) = 0, (4) 


where U denotes either p, pe, or put, and F, and F„ denote the inviscid and viscous contributions to 
flux. Equation (4) is transformed to a local computational coordinate system for each discrete element. 
(£,r],C) denote the local coordinate system with Jacobian J = d(x±, X 2 , xz)/d(£, 77, £). Equation (4) has 
same form in the local coordinates: 


dU 

~dt 


+ V ■ (Fi 


F„) = 0, 


the 

Let 

the 

( 5 ) 


but with U = U\J\, Fi = J-'Fil J\, and F v = J- 1 F„|J|. 
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A. DG Discretization 


The implementation of DG used in this study follows the quadrature-free form described in Refs. 4,6,7, and 
is described only briefly here. The computational domain is subdivided into nonoverlapping elements that 
cover the domain. The DG discretization is formulated locally in each element in a similar manner. The 
solution within each element is approximated as an expansion in a local basis set {b n }, usually polynomials 
of degree < p, 

N p 

U — ^ ' b n u n , 

n — 0 

where N p denotes the number of terms in the basis set of degree p. The current work uses three-dimensional 
monomials of the form for all ( i,j,k ) triplets satisfying 0 < i + j + k < p. However, both the 

formulation and the simulation code readily apply to lower dimensional formulations simply by constraining 
the basis to lower dimensions. The number of unknowns in each element is the number of physical variables 
times the size of the basis set, N p . An equal number of equations governing these unknowns is derived by 
multiplying the governing equations by each member of the basis set, and integrating over the element. The 
integrals of the flux terms are integrated by parts to obtain the following weak form: 

[ bj^-dn- I • (Fi - F„) dn + [ b n {Fi - F„) • n ds = 0 Vfc : 0 < k < N p , (6) 
Jn vt J n J dn 

where denotes the element, dQ denotes the element boundary, and n denotes the outward unit normal. 
Following the quadrature-free formulation, the fluxes are expanded in a similar manner: 

M M 

Fi = ^2 bnfi ’ k and F „ = ^2 b n fv,k- 

fc= o fc= o 


The degree of the flux expansions is allowed to be higher than that of the solution with the result that 
M > N p . The element boundary integral is evaluated piecewise on discrete segments or edges, dClj, where 
each edge is shared with a neighboring element. The two neighboring elements on either side of an edge 
segment have independent local approximations for the solution and fluxes. To resolve this ambiguity, the 
fluxes in the edge integrals are replaced by numerical fluxes that are a function of the solutions in both 
neighboring elements. The numerical edge fluxes are represented in terms of a lower dimensional edge basis, 
bk,j- 

M M 

Fi ‘ n la! 2 J - — k Fj = ^ ' bkjfj,k,j and F„ • n| af! ^ — > F v = ^2 bk,j fv,k,jt 
k = o fc= o 


where dttj denotes an individual edge and M denotes the number of terms in the edge basis. 

The inviscid edge flux is modeled by an approximate Riemann flux. The local Lax-Friedrichs flux is 
simple and inexpensive to implement, and has worked well with DG for smooth flows. 4, '~ 9 However, the 
HLL flux 10 has been shown 5 to be more robust accurate for flows with shock waves, and is used in all 
simulations presented here. 

The solution gradients required for the viscous terms are evaluated using the DG methodology in a similar 
manner. Let 

N p 


& w 



= Vw, 


k=0 


where w denotes either Ui or T. It immediately follows that 


— / S7b n w dfl + / b n w n ds = 0. 
Jq J dQ 


As before, the multi-valued edge flux is replaced by a numerical flux, icn| sri = {on. Both F v and foil are 
evaluated as described in Ref. 4. 

Because all of the fluxes are represented as expansions in the basis set, the integrations can be performed 
directly to produce a matrix equation of the form: 


M 


du n 

dt 


- v • [fi.*, — f^/s] + '22 Bj 


^ i,k,j ^v,k,j 


= 0, 
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where M, V, and B are derived in Ref. 6 along with further details of the quadrature-free formulation. 

Reference 5 investigated a number of implementation details and variants of DG that influence the 
accuracy and robustness of the method for flows with shock waves. In most cases, these variants improve 
the accuracy of DG when shocks are not present, However, the exception is the strategy of augmenting the 
viscosity in the neighborhood of the shock wave to improve the robustness of the method. While the earlier 
work demonstrated that the technique did indeed improve robustness, it also revealed some undesirable side 
effects. In particular, for the preliminary axisymmetric simulations of the current flow, the forcing could be 
increased by an order of magnitude over what was possible with out applying the technique. However, shock 
sensors typically used by the approach would also be triggered by the disturbance itself, locally increasing 
the viscosity and potentially corrupting the simulation. Due to the developmental nature of the augmented 
viscosity technique, it was decided not use the technique in the 3D simulations presented here. While this 
choice did restrict the maximum forcing amplitude that could be run, the more gentle nature of the long- 
timescale forcing (Eq. 9 below) allowed 3D simulations with sufficiently large amplitudes to run successfully. 

III. Simulation 

Time accurate unsteady simulations were performed starting from a shock- free steady flow solution 
described in Ref. 11. The initial residuals were subtracted from all following values to eliminate small 
transients that could arise due to differences in the discretization schemes and grid resolution. The Reynolds 
number is 9.55 x 10 5 6 / m. The discretization is a fourth-order DG on tetrahedral elements based on the 
quadrature- free formulation. 6, 12 Time advancement uses the explicit third-order TVD-RK scheme of Shu. 13 

A. Computational Domain 

The computational domain, shown in Fig. 2, is defined by the tunnel wall, two nearly conical sections, two 
ramp sections, and two x=constant planes for the inflow and outflow boundaries. The conical and ramp 




(a) Domain boundaries (b) Grid cross-section on tunnel wall and several cross-flow 

sections; showing every fourth line. 


Figure 2. Computational domain and grid. 


sections are inclined steeper than the local Mach cones expected in the Mach 6 freestream such that these 
inflow and outflow boundaries are not of a mixed type (i.e. characteristic directions at boundary points are 
either all inbound or all outbound). The underlying grid is a curvilinear, single-block structured grid in 
which each quad-element is subdivided into five tetrahedrons. Grid stretching is applied in the wall-normal 
direction to place approximately half of the points within the boundary layer. Within the boundary layer, 
grid points are clustered toward the forcing region and toward the centerline. In the streamwise direction, 
the domain extends from x = 1.85m to x = 2.584m, for most of the results show here. The downstream 
extent of the domain in the current simulations is limited by the domain of the steady starting flow solution 
taken from Ref. 11. A typical grid size of 345 x 54 x 32, with five tetrahedrons per quad and 20 degrees 
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of freedom per polynomial of degree three yields approximately 60 million degrees of freedom per equation. 
The accuracy of simulation results was evaluated by comparing them to results on a coarser grid for isolated 
cases, and is judged to be more then adequate for the current flow. These comparisons will be shown in a 
later section. 

B. Modeling Effect of Electrode 

The primary effect of the electrode is to inject energy into the flow field. This effect is modeled numerically 
by adding a source term of the form y, z)T(t ) to the energy equation (right hand side of Eq. 2). Little 

was know at the outset of this work about the forcing produced by the experimental apparatus; and in the 
preliminary two-dimensional axisymmetric simulations, 5 the source term was modeled as simple pulse: 

S(x, y, z) = [(1 + cos( 7rr))/2] 2 (7) 

T{t) = sin( 7rt ) 4 (8) 


where 

r = min(l, y/(x- x s ) 2 + {{y - y s )/a) 2 + ((z - z s )//3) 2 / r f ), 
t = min(l, t/td), 

td is the duration of the forcing, f is a directionally scaled distance from the center of the source, rf is the 
radius of the forcing region, and Aq is the forcing amplitude. The origin of the forcing is at (x, y , = 

(1.924, 0.117, 0.0). The radius of the forcing region was set to 0.00127m, which approximately equals the 
electrode gap. To elongate the forcing region in the streamwise direction, a and /? are set to 1/2. The 
duration was set to 2 x 10 _5 s by matching the period of the disturbance at x = 2.201. The amplitude Aq 
was tuned to match the disturbance amplitude at multiple probe locations and typically ranged between 
5 x 10 4 to 6 x 10 5 . 

Recent experimental work 3 suggests that the duration of the forcing is much longer than originally 
assumed, and the following model for T(t) is used in the 3D simulations. 

T(t) = 16 [sin( 7 rf ) + sin( 2ni ) 2 /2] [(1 — Ai) + A x cos( 7r( 2 1 — td/3 )/ t ^)] / 27. (9) 

The primary forcing is now asymmetric, and the new parameters, A\ and th denote the amplitude and 
period, respectively, of a short-period component that is matched to the period of the disturbance. Forcing 
parameters td = 1.624 x 10 -4 (s), th = 2.17 x 10 _5 (s) and rf = 0.01m are used in the simulations presented 
below; A® and A\ are varied to approximately match the character of experimental results at x = 2.055 and 
the period of the disturbance at x = 2.201. 

C. Axisymmetric Simulations 

Axisymmetric simulations were used to both size the forcing parameters and to evaluate mesh size require- 
ments. Figures 3(a) and (b) illustrate the two types of forcing functions and the disturbances they produce. 
Pressure data at each probe are normalized by their local initial (undisturbed) values, Pq- The simple 
sin(x) 4 type of forcing produces a symmetric wave packet at a targeted frequency. However, the experimen- 
tal measurements indicated that, at this location, the disturbance should be small and superimposed onto a 
long-period transient. The more complex asymmetric forcing given by Eq. 9 has this desired characteristic. 
The amplitude of the long-time transient seen in the experiment is still much larger than that of the simu- 
lation. However, further effort to reverse engineer the simulation forcing is beyond the scope of this work. 

Figures 4(a) and (b) show the results at two probe locations for a range of mesh sizes. Results are grid 
converged on the 320 x 75 grid, and are acceptable on the 213 x 50 mesh. Most 3D simulations are performed 
on a grid similar to the 320 x 75 2D grid. However, the number of points in the wall normal direction are 
decreased in the 3D case by increasing the stretching outside the boundary layer, and by decreasing the 
distance to the outer boundary. The number of points in the stream wise direction was increased in the 3D 
grid, relative to the 2D grid, to provide better resolution near the forcing. 
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(a) Temporal forcing function T(t): blue line indicates the sim- 
ple symmetric pulse, black lines indicated the complex asym- 
metric pulse given by Eq. 9 for given values of A\. 

Figure 3. Forcing functions 



(b) Disturbance produced by the two types of forcing at x = 

2 . 201 . 


and disturbances produced. 




Figure 4. Mesh refinement study at two probe locations. 
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D. Three-Dimensional Simulation Results 


The experimental data consist of the time history of pressure measured by probes on the tunnel wall at 
several locations downstream of the pulse glow perturber. The data are post-processed as described in 
Ref. 1. Only 5 of the probe locations fall within the computational domain of the simulation. Similar data 
were collected during the simulations at each of these locations, as well as at an additional 478 points. The 
forcing parameters, Aq and A\ were tuned to approximately match the disturbance amplitude observed in 
the experiment at x = 2.201m. Figures 5 (a) and (b) show the experimental and simulation results for the 
common set of probe locations with A$ = 3,710 and A\ = 0.078167. Qualitatively, the principal period of 
the disturbance and its group velocity agree well with the experiment. In the experiment, the wave packet 
has a short period component superimposed over a much broader transient. The simulated results have a 
similar feature, but it is much less pronounced. 



t(s) 


(a) Experimental results. 


0.3 - 

x=2.055 m 

x=2.201 m 

x=2.302 m 

0.2 - x=2.378 m 


x =2.480 m 



_n 0 I 1 i : i : : : ! 1 1 1 1 1 1 . I 

0.0002 0.0004 0.0006 0.0008 0.001 

t(s) 

(b) Simulation results with Aq = 3,710, A\ = 0.078. 


Figure 5. Time history of pressure at probe locations common to both data sets. 


To precisely distinguish between the two components, the data are further decomposed into a long- 
timescale transient and a short-timescale symmetric disturbance. Figures 6 (a) and (b) illustrate the results 
of this process for the experimental data at the probe location x = 2.480. The decomposition method is a 
heuristic procedure based on identifying extrema. The method works well when the mode scales are well 
separated but can locally fail under some well understood circumstances. For the current purposes, the 
method can be reliably used to precisely quantify features of the symmetric disturbance such as its period 
and growth rate. The decomposition is constructed by creating envelopes defined by cubic fits through 
local maximums and minimums. The long-timescale transient is defined as the average of the upper and 
lower envelopes. The short-timescale symmetric disturbance is defined as the difference between the actual 
signal and the long-timescale component. Subtracting the long-timescale from the upper and lower envelopes 
yields the envelope of the short-timescale symmetric disturbance. In the following discussion, the disturbance 
amplitude A d is defined to be the maximum of the disturbance envelope, £p ca k is the time at which the 
disturbance envelope reaches its peak, and the principal period ^periocl the peak-to-peak time interval 
averaged over the three cycles nearest the peak disturbance. Growth rates per meter, g : are computed from 
adjacent probes located at Xi and Xi-i as: 

g = \a(A di i/A d 'i-i)/{xi - xi- 1 ) 


Applying this procedure to all data samples gives the results shown in Fig. 7. The amplitude of the 
long-timescale transient is larger than the disturbance for the first three probe locations of the experimental 
results. At the most upstream probe location, the disturbance is nearly undetectable in the full pressure 
trace, and the decomposition method is not reliable for this probe. In the simulation results, the long-scale 
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Figure 6. Decomposition of data into long-timescale transient and short-timescale symmetric disturbance. 


transient amplitude is comparable to that of the disturbance only at the most upstream probe, and the 
decomposition procedure is judged to be reliable at all probe locations. 




(a) Experimental results. (b) Simulation results with Ao = 3,710, A i = 0.078. 

Figure 7. Time history of long-timescale transient and short-timescale symmetric disturbance at probe loca- 
tions that are common to both sets. 

With the disturbance component isolated, disturbance features such as amplitude, growth rate, group 
velocity, and period can be precisely determined, and are presented in Figs. 8 and 9. The simulation 
data include additional probe locations between those of the experiment; and the experimental data include 
additional probe locations that are downstream of the simulation domain. The figures show results for several 
combinations of forcing amplitude. However the simulation results agree very well with the experimental 
result only for the case with the largest forcing amplitude (black curve). 

It was expected at the outset that for any given value of long-mode forcing A 0: the short-mode forcing 
parameter A\ could be tuned to match the experimental results; however, such was not the case. Fig- 
ures 8 and 9 also show results from two additional simulations in which the forcing amplitude of the short 
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(b) Growth rate (per meter) of disturbance. 

Figure 8. Amplitude and growth rate of the symmetric disturbance. 




Figure 9. Time of peak disturbance amplitude and period of disturbance at peak. 
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mode (AqAi) is held fixed, while the amplitude of the long mode (A)(l — Ai)) is reduced by half (blue 
curve), and half again (green curve). In these cases with lower long- mode forcing, the disturbance amplitude 
is matched at x = 2.055, but the growth rate is nearly constant through the computational domain, and is 
lower than the peak values observed in both the experiment and in the simulation with higher long-mode 
forcing. Figure 9(a), which plots the time at which the disturbance amplitude reaches its value peak, tp ea k, 
indicates that the speed of the disturbance is not significantly influenced by the level of the long-mode 
forcing. All simulations give similar values for the disturbance period. However, the simulation with the 
largest long- mode forcing does shift more towards the experimental values (e.g., lower upstream values that 
gradually increase in the streamwise direction). 


IV. Linear Stability Analysis 

To gain insight into the differences seen in the simulations with the variations the forcing amplitude, 
a linear stability analysis was performed using the Langley Stability and Transition Code (LASTRAC). 14 
The mean flow input required by LASTRAC was constructed from the same mean flow used in the flow 
simulations. Fig. 10(a) shows growth rates predicted by LASTRAC for an analysis of second-mode instability 
wave over a narrow range of frequencies. Each mode has a peak growth rate of about eight, which is much 
lower than the rates of the experiment, or of the simulation with large forcing amplitude. Also, none of the 
modes have a sustained growth rate over the entire spatial range. This suggest, that the forcing in both the 
experiment and in the simulation must create a disturbance with a rich complement of modes. Integrating 
the growth rates of each mode gives the disturbance amplitudes shown in Fig. 10(b). The thick dashed line 
is the result of averaging the amplitudes of all modes that are growing in the region of interest. The growth 
rate of the average is about four, shown as the thick dash line in Fig. 10a and agrees well with the simulation 
results with low amplitude forcing. 


20 - 
g(1/m) ; 
15 - 

10 - 


f = 52,500Hz 
f = 50,000Hz 
f = 47, 500Hz 
f = 45,000Hz 
f = 42, 500Hz 
f = 40,000Hz 
f = 37, 500Hz 
f = 35,000Hz 
Average 



2.2 2.3 

X(m) 

(a) Growth rates. 



2.2 2.3 

X(m) 

(b) Integrated amplitudes. 


Figure 10. Results from linear stability analysis on a uniform spectrum. 

It is possible to expand the numerical forcing function to reveal all of the modes it produces. Expanding 


11 of 13 


American Institute of Aeronautics and Astronautics 


a simplified version of trigonometric functions of Eq. 9 gives: 

[sin( irt/td ) +| sin( 2nt/t d ) 2 ] [(1 - A±) + A 1 cos( 2nt/t h )] 

= (1 — A{) {5 + 4 [cos( nt/t d ) — cos( 2nt/t d ) — cos( 3nt/t d )] — cos( 47 rt/t d )} / 8 
+A 1 { 10cos( 2 7 rt/th ) — cos( 2 Trt(2/t d — 1 /th) ) — cos( 2 Trt(2/t d + 1 /th) ) 

+ 4 [ COs( TTt(l/t d - 2/t h ) ) + COs( TTt{l/t d + 2/th) ) 

- cos( nt(3/t d - 2/th) ) + cos( nt(3/t d + 2/th) ) 

- cos( 7rt(4/ t d - 2/t/j) ) + cos( irt(A/t d + 2/t/,) )]}/16 

(10) 

Evaluating the above using t d = 1.62451 x 10 4 s and th = 2.171158 x 10 5 s reveals that the forcing is exciting 
the following frequencies: 3,078Hz, 9,233Hz, 6,155Hz, 12,311Hz, 33,745Hz, 36,825Hz, 39,903Hz, 42,980Hz, 
46,058Hz, 49,136Hz, 52,214Hz, 55,292Hz and 58,370Hz. Figure 11 shows the growth rates and integrated 
amplitudes predicted by LASTRAC for the unstable modes of this list. The thick black lines indicate the 
results from uniformly averaging the individual modes and of a weighted average using the actual weights 
suggested by Eq. 10. 
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(a) Growth rates. (b) Integrated amplitudes. 

Figure 11. Results from linear stability analysis on the spectrum produced by the simulation forcing. 


V. Conclusions 

Simulations of induced disturbances in a laminar Mach 6 boundary layer have been performed using a 
high-order DG method on an unstructured grid. At a high forcing amplitude, the amplitude, growth rate, 
group velocity and period of the disturbance agree well with experiments performed in the Boeing/AFOSR 
Mach-6 Quiet Tunnel at Purdue University. However, the peak growth rates of both cases are three times 
greater those predicted by the linear stability analysis. 

In simulations with low forcing amplitude, the growth rate of the disturbance is lower and nearly constant 
over the domain. However, the group velocity and period of the disturbance are still similar to that of 
the experiment and of the simulation with high forcing amplitude. Linear stability analysis predicts that 
frequencies between 35kHz to 55kHz are unstable over some portion of the domain of interest, but that no 
single mode is unstable over the entire region. Individual modes have peak growth rates of about eight. Both 
uniform and weighted sums of unstable modes produce growth rates of about four throughout the region, 
which agrees well with simulations using lower forcing amplitudes. 

The high growth rates seen in the experiment and in the simulations with large amplitude forcing are 
most likely due to nonlinear flow effects. The uneven growth rate seen in the simulation is attributed to the 
forcing function that populates only a few discrete modes, where as the pulse glow perturber likely induces 
a continuous spectrum of modes. The sensitivity of the growth rate to the overall forcing amplitude, even 
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while holding the amplitude of the short-period forcing mode fixed, suggest that details of the forcing play 
a critical role in the outcome. Additional measurements to characterize the pulse glow perturber and its 
effects on the immediate flow field are needed for rigorous CFD validation. 
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